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We investigated the nematic to smectic transition undergone by parallel hard spherocylinders in 
the framework provided by the residual multi-particle entropy (RMPE) formalism. The RMPE is 
defined as the sum of all contributions to the configurational entropy of the fluid which arise from 
density correlations involving more than two particles. The vanishing of the RMPE signals the 
structural changes which take place in the system for increasing pressures. Monte Carlo simulations 
carried out for parallel hard spherocylinders show that such a one-phase ordering criterion accu- 
rately predicts also the nematic-smectic transition threshold notwithstanding the almost continuous 
character of the transition. A similar quantitative correspondence had been already noted in the 
case of an isotropic fluid of freely rotating hard spherocylinders undergoing a transition to a nematic, 
smectic or solid phase. The present analysis confirms the flexibility of the RMPE approach as a 
practical and reliable tool for detecting the formation of mesophases in model liquid- crystal systems. 

PACS numbers: 64.60.-1, 64.70. Md, 61.20.Ja, 65.40.Gr 



I. INTRODUCTION 

Recently we investigated the formation of liquid- 
crystalline mesophases in a fluid of hard spherocylinders, 
i.e., cylindrical segments of length L and diameter D 
capped at each end by a hemisphere of the same diame- 
ter 0] . More specifically, we analyzed the ordering of the 
homogeneous and isotropic fluid into a nematic, smec- 
tic or solid phase within the framework provided by the 
one-phase entropy-based criterion originally proposed by 
Giaquinta and Giunta 0- This criterion can be imple- 
mented through the well known multi-particle correlation 
expansion of the configurational entropy i^J : 

oo 
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In the above formula Sox is the excess entropy per particle 
in units of the Boltzmann constant and Sn is the n-body 
entropy that is obtained upon re-summing spatial corre- 
lations between up to n particles. In particular, the pair 
entropy per particle of a homogeneous and isotropic fluid 
of nonspherical molecules can be written as Q: 

- J{g{r,u;')\n[g{r,u;')] 

-g{v,u^) + l}dvdLu\ (2) 

where g{r,Lu'^) is the pair distribution function (PDF) 
which depends on the relative separation r between two 
molecules and on the set of Euler angles = [ti;i,u;2] 
that specifies their absolute orientations in the labora- 
tory reference frame. The quantity Q represents the in- 
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tegral over the Euler angles of one molecule, while p is 
the particle number density. 

In order to identify the ordering threshold of the fluid, 
we monitor the behavior — as a function of the number 
density — of the so-called residual multiparticle entropy 
(RMPE). This quantity is defined as the difference 

As = Sox - 32 ■ (3) 

At variance with the pair entropy, the RMPE exhibits a 
nonmonotonic behavior as a function of p. In particular, 
it is negative at low densities, and becomes positive as 
a more ordered phase is approached. The relevance of 
the condition As = as a one-phase ordering criterion 
has been documented for a variety of phase transitions, 
both in continuous fluids as well as in lattice-gas model 
systems 0. 

Even for systems composed of non-spherical molecules, 
such as hard spherocylinders, it turns out that the 
ordering thresholds detected through the zero-RMPE 
condition systematically correlate with the correspond- 
ing phase-transition points, whatever the nature of 
the higher-density phase coexisting with the isotropic 
fluid Ul ■ Cuetos and coworkers have successfully applied 
the RMPE criterion to hard spherocylinders with an at- 
tractive square well and to spherocylinders with a soft 
repulsive core 0|. 

In this paper we intend to analyze the predictions of 
the RMPE approach in the case of the nematic-smectic 
transition undergone by parallel hard spherocylinders 
with aspect ratio L/D = 5. The onset of smectic order 
out of a nematic phase represents the next step in the 
process which, upon compression of the isotropic fluid, 
eventually leads to the formation of the fully crystalline 
solid. The phase behavior of this model has been investi- 
gated with several numerical simulations B 13 as well 
as theoretical studies (see e.g. 0,0, 0| and references 
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contained therein) over the whole L/D range. Such stud- 
ies have ascertained the existence of a second-order phase 
transition from a low-density nematic state to an inter- 
mediate smcctic phase. On the other hand, the general 
features of the higher-density region of the phase diagram 
has long been debated as far as the stability a columnar 
phase is concerned 0. 

The paper is organized as follows: In Sections 2 and 3 
we describe the theoretical framework and the numerical 
simulation technique. The results are presented in Sec- 
tion 4 while Section 5 is devoted to concluding remarks. 



II. THEORY 

The formalism developed in Q for the isotropic fluid 
needs to be modified in order to describe the nematic and 
smectic mesophases. Even for uniaxial molecules, the 
numerical computation of the full pair distribution func- 
tion g{Y,Ljp) is a formidable task. Costa and coworkers 
have already shown that, in order to reproduce the phase 
boundaries of the isotropic fluid, it is enough to take into 
account the dependence oig{Y, lo^) on the centers-of-mass 
separation r and on the angle 9 formed by the molecular 
axes of two spherocylinders This angle shows up as 
the critical parameter which accounts, by itself, for the 
reduction of orientational states that is ultimately offset 
by the gain of translational entropy at high densities |l3j| . 
On the other hand, it is rather obvious that the angle 9 
conveys no useful information on the structural process 
that may eventually lead to the formation of a smectic 
phase out of a nematic phase. 

In order to resolve ordering effects associated with the 
orientational degrees of freedom on one side, and with the 
modulation of the density along the nematic director on 
the other side, we shall restrict our analysis to a system 
composed of parallel spherocylinders. Assuming that the 
molecules are perfectly aligned does not significantly alter 
the main aspects of the phenomenology that we want to 
investigate. As emphasized above, different mechanisms 
drive the phase transformation of an isotropic or of a ne- 
matic fluid into a more ordered phase. Moreover, at the 
nematic-smectic transition threshold the nematic fluid is 
typically characterized by a very high degree of relative 
alignment of the molecules. For example, the nematic 
order parameter is about 0.85-0.89 for spherocylinders 
with aspect ratio L/D — 5 14] . The constraint on the 
relative orientation of the particles hugely simplifies the 
expression of the PDF. In fact, in order to specify the 
spatial configuration of a given spherocylinder relative 
to that fixed at the origin one just needs two parame- 
ters: the distance r and the angle d formed by r and by 
the nematic director. The other polar angle (f> is actually 
averaged out on account of the macroscopic cylindrical 
symmetry of the model. We emphasize that in this case, 
at variance with the freely rotating model already inves- 
tigated in 0, no approximation is needed in order to 
compute the PDF of the system. Correspondingly, the 



pair entropy of the fluid can be written as: 
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{g(r,z9)ln[g(r,,9)]-g(r,i9) + l}. (4) 

Using a formalism analogous to that introduced in jj], 
we can extract from Eq. the excluded-volume contri- 
bution to the pair entropy of the fluid. This contribution 
arises from the space integration carried out over the re- 
gions where g{r, d) — 0: 



S2 = -B2P + s'2 



(5) 



In Eq. ||SJ| B2 is the second virial coefficient of hard par- 
allel spherocylinders that is just four times the volume, 
V\isci of a spherocylinder: 



B2 ^ 4[iTr/A)D^L + (7r/6)L»^ 



(6) 



The residual contribution, S2, explicitly accounts for the 
decrease of the pair entropy associated with the onset 
of interparticle correlations at short and medium-range 
distances. For parallel spherocylinders, a further separa- 
tion of this term into a translational and an orientational 
contribution as done in 1] does not make much sense. 

The total excess entropy of the nematic fluid can be 
evaluated upon integrating the equation of state (EoS): 
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In Eq. Q P is the pressure, [3 is the inverse temper- 
ature in units of the Boltzmann constant, and p is the 
density of a suitably chosen reference state. The EoS of 
the fluid was obtained through a Monte Carlo sampling 
of the system at different pressures, while the quantity 
5cx(p) was calculated using the Widom test-particle in- 
sertion method [isl IT^ in order to calculate the excess 
chemical potential, /Xex, at low enough densities. The 
excess entropy can then be obtained through the ther- 
modynamic equation: 



= -PMcx H 1 

P 



III. SIMULATION 



(8) 



We investigated the phase diagram of a system of par- 
allel hard spherocylinders with elongation L/D = 5, 
spanning a density range which goes from a rather dilute 
nematic state up to the smectic transition threshold. We 
carried out Monte Carlo (MC) simulations at constant 
pressure as is usual for systems of non-spherical hard-core 
particles where it may be difficult to calculate the equa- 
tion of state in a constant-volume simulation through the 
contact values of the distribution functions [l^. 

The typical sample was composed of iV = 1500 par- 
ticles, aligned along the z axis and enclosed in an or- 
thorhombic box with edges L^ = Ly = l/SL^. In order 
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to quantify the influence of the size of the system as well 
as of the shape of the simulation box, we also investi- 
gated the behavior of a system composed of 500 particles 
aligned along the main diagonal of a cubic box and of a 
fluid of 768 particles enclosed in a orthorhombic cell with 
the same relative dimensions used for the main sample. 

All thermodynamic states at constant pressure were 
sequentially generated from a translationally disordered 
low-density configuration upon gradually compressing 
the nematic fluid. The equilibration period was typi- 
cally 10^ MC cycles, a cycle consisting of an attempt 
to change sequentially the center-of-mass coordinates of 
each molecule followed by an attempt to modify the vol- 
ume of the sample. Simulation data were obtained by 
generating chains consisting of 5 x 10^ — 20 x 10^ MC 
cycles, depending on the pressure. Equilibrium aver- 
ages and standard deviations were computed by dividing 
chains into independent blocks. During the production 
runs, we cumulated different histograms of the PDF. In 
particular, g{r, •§) was sampled at intervals Ar and Ai9 of 
0.05D and 1°, respectively. Different choices of Ar and 
At? were investigated for P* — 1.0 and P* = 2.0. As 
for the Widom insertion method, 100 trial insertions per 
MC cycle turned out to be sufficient to insure a stable 
statistics for the excess chemical potential of the fluid. 

In the presentation of the results, we shall refer to 
the reduced density p* = p/pcp, where Pcp = 2/[\/2-l- 
(L/D)V3] is the maximum density attained by paral- 
lel spherocylinders at close packing, and to the reduced 
pressure P* = pPV^sc- 



IV. RESULTS 

The current MC results for the EoS of the model are 
presented in Fig. ^ The comparison between the data 
collected for 1500 particles (in an orthorombic box) and 
those for 500 particles (in a cubic box) shows that the 
EoS of the fluid is rather weakly affected by the size of 
the sample only at high densities. The results for 768 
particles in an orthorombic box were not reported in 
the graph since they can be hardly resolved from those 
pertaining to the largest size investigated. The present 
EoS is in good agreement with that computed by other 
authors Q over the whole density range explored, 
apart from a modest deviation, at high densities, from 
the molecular dynamics results obtained by Veerman and 
Frenkel . Both these authors as well as Stroobants and 
coworkers noted a weak change of slope in the EoS at 
the transition that was estimated to occur for p* = 0.46. 
Other structural properties of the system were evaluated 
in such as the distribution functions associated with 
parallel and perpendicular correlations, below and above 
the observed cusp. Given the absence of either transla- 
tional order within the layers or hysteresis effects asso- 
ciated with a compression/expansion cycle, Stroobants 
and coworkers concluded that the nematic phase actu- 
ally transforms into a smectic phase through a continuous 



phase transition. Veerman and Frenkel later estimated 
the transition density by studying the critical slowing 
down of the intermediate scattering function evaluated 
along the nematic director Their result agrees with 
that reported by Stroobants and coworkers f^. 

We verifled that, during the simulation, the system 
spontaneously transforms into a smectic phase, as can be 
seen through the snapshots of equilibrated configurations 
that are shown in Fig. [3 Such snapshots were taken for 
densities across the transition threshold. Traces of smec- 
tic order can be hardly detected for P* — 2.00, a pressure 
corresponding to an average density equal to 0.44. The 
ordering of the fluid is already apparent just beyond the 
transition threshold for P* = 2.30 and p* ~ 0.47. Fully 
developed arrangements of several well-separated smectic 
layers distributed along the z direction are clearly visible 
for P* ~ 2.50 and p* ~ 0.49. Instead, the transversal 
arrangements of the spherocylinders along the x-y plane 
look disordered in all states investigated (compare left 
and right panels in Fig. [21). 

We show in Fig. O the excess entropy of the fluid plot- 
ted as a function of the reduced density. The data were 
obtained using both the thermodynamic integration of 
Eq. Q as well as Widom's ghost-particle method imple- 
mented at moderately low densities. We also show the 
datum recently reported by Koda and Ikeda joj which 
was obtained using a multistage Widom test, based on 
the gradual insertion of a ghost particle with a variable 
shape. The modest discrepancy between this finding and 
the present results is likely due to the multistage method, 
that is known to work better at high densities. On the 
other hand, our results for the excess entropy are closely 
interpolated by the expression obtained upon integrating 
the five-term virial fit of the EoS reported in 0. 

The PDF of the fluid, g{r, d), was plotted in Fig. Has 
a function of -d for a set of interparticle distances at in- 
creasing pressures across the transition point [P* ~ 2.2). 
We note, for separations r < L + D, the existence of a 
forbidden range around iS — and tt that is due to the 
overlap of two spherocylinders. This correlation gap de- 
creases with increasing intermolecular separations. The 
maximum attained by g(r, i?) corresponds to the hard- 
core contact between spherocylinders and its height con- 
sistently decreases with increasing interparticle distances. 
For r > L + D, the entire angular range between and 
TT can be sampled by a second spherocylinder. Of course, 
there is a blow up of the contact value of the PDF for 
"i? = 0, TT, which progressively decreases for increasing dis- 
tances. An increase of the pressure enhances the overall 
structure of the PDF; however, no speciflc signature of 
the nematic-smectic transition can be detected in the re- 
sulting modiflcation of the PDF. 

Upon plugging the PDF into Eq. Q), we obtained the 
pair entropy that was plotted together with the excess en- 
tropy in Fig.|5las a function of the density. The RMPE 
exhibits a change of sign from negative to positive values 
for p* ~ 0.453. This threshold practically coincides with 
the currently accepted nematic-smectic transition density 
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{p* = 0.46) HHl. The simulations carried out for 500 
particles enclosed in a cubic box show a similar behav- 
ior of the RMPE whose morphology is not substantially 
modified with respect to that of 1500 spherocylinders in 
an orthorombic box. This finding is noteworthy in that 
it demonstrates the high sensitivity of the RMPE to the 
structural changes occurring in the fluid, irrespectively 
of the difficulty to accommodate a well resolved smectic 
layering in a relatively small cubic box. 

Figure |H| shows the quantity that, upon integration, 
yields the pair entropy in Eq. I^J. For very short inter- 
particle distances the "i? range which can be sampled by 
two neighboring particles is rather limited. As a result, 
angular correlations are strong and give rise to the deep 
well in the integrand function. A steady modulated in- 
crease follows up to r = L -t- 13 where a negative jump 
witnesses the onset of new strong correlations between 
spherocylinders lying on top or below the central refer- 
ence particle. This effect is the integrated counterpart 
of the behavior that was already discussed for g(r, ??) in 
Fig.EI 

The pair entropy was finally resolved into an excluded- 
volume and a correlation term in Fig. [7| (see Eq. |3J| . It 
is clear that the second-order virial term cannot account 
by itself for the crossover between Scx and S2 that is in- 
dicative of the smectic ordering of the nematic fluid. 

As for the numerical reliability of the current results, 
we note that the statistics cumulated on g(r, ??) was 
such that a smooth integration over the sampled points 
yielded stable values for S2, the dispersion being always 
lower than a percent of the average value. However, we 
observed a moderate sensitivity of the angle-dependent 



quantities in Eq. on the resolution of the angular 
width A-ff. In order to gain a better insight into the 
numerical accuracy of the calculations, we performed a 
series of test runs for two different pressures (P* = 1.0 
and 2.0) with several r and i? grid meshes. The tests 
were reported in Fig. |S1 where the higher sensitivity of 
S2 on the angular grid size — as compared with the ra- 
dial one — is quite manifest. We estimated, in the limit 
of Ar, Ad ^ 0, an asymptotic shift of S2 toward lower 
values not larger than 2%. As a result, we expect a com- 
parable shift of the transition threshold signalled by the 
RMPE to lower densities. 



V. CONCLUDING REMARKS 

In this paper we have analyzed the residual multipar- 
ticle entropy (RMPE) of parallel hard spherocylinders 
with aspect ratio L/D = 5 across the nematic-smectic 
transition. The correspondence between the intrinsic or- 
dering threshold detected through the vanishing of the 
RMPE (p* ~ 0.453) and the independently ascertained 
phase-transition density (p* — 0.46) is quantitative. It 
also turned out that the indication of the zero-RMPE cri- 
terion is not significantly affected by the size and shape 
of the simulation box, even when the smectic layering of 
the fluid cannot be easily accommodated in the sample 
as is the case of a small cubic box. This finding further 
corroborates the sensitivity of the RMPE to the ordering 
of the fluid on a local scale. 
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FIG. 1: Equation of state of parallel hard spherocylinders with aspect ratio L/D = 5. Line with solid circles and crosses: 
this work for A'^ = 1500 and A'' = 500 molecules, respectively. Triangles, squares and diamond: simulations by Veerman and 
Frenkel |3, Stroobants and coworkers 0, and Koda and Ikeda Q, respectively. The error bars are systematically smaller than 
the size of the markers. 




FIG. 2: Snapshots of equilibrated configurations for (from top to bottom) P* = 2.50 ((p* = 0.23)), P* = 2.30 ({p* = 0.44)), 
P* = 2.00 ((p* = 0.47)), and P* = 0.50 ((p* = 0.49)). Left panels: projections of the centers of mass of the spherocylinders 
onto the x-y plane. Right panels: configurations along the z axis. 
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FIG. 3: Excess entropy plotted as a function of the reduced density. Line with solid circles: this work with = 1500; diamonds: 
Widom insertion-method estimates. Square: multistage Widom test 0. Dashed line: excess entropy evaluated by integration 
of the five-term virial fit of the simulation data of Ref.Q. The vertical line indicates the nematic-smectic transition threshold 
according to Ref. and |3- 
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FIG. 4: Pair distribution function g{r,'d) plotted as a function of '& for different distances r/D and for increasing pressures. 
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FIG. 5: Residual multiparticle entropy (circles) resolved into the excess (squares) and pair (solid diamonds) contributions. 
Lines are smooth interpolations of the simulation data. 




FIG. 6: Integrand function appearing in Eq. @, after integration over plotted for several pressures. 
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FIG. 8: The pair entropy S2 plotted as a function of the grid mesh parameter Ai? and for several choices of Ar, for P* = 1.00 
(solid symbols) and P* = 2.00 (open symbols). Triangles, Ar — O.IOD; circles, Ar = 0.05-D; squares Ar = 0.02_D. 



